Time‐calibrated phylogeny and ecological niche models indicate Pliocene aridification drove intraspecific diversification of brushtail possums in Australia

Abstract Major aridification events in Australia during the Pliocene may have had significant impact on the distribution and structure of widespread species. To explore the potential impact of Pliocene and Pleistocene climate oscillations, we estimated the timing of population fragmentation and past connectivity of the currently isolated but morphologically similar subspecies of the widespread brushtail possum (Trichosurus vulpecula). We use ecological niche modeling (ENM) with the current fragmented distribution of brushtail possums to estimate the environmental envelope of this marsupial. We projected the ENM on models of past climatic conditions in Australia to infer the potential distribution of brushtail possums over 6 million years. D‐loop haplotypes were used to describe population structure. From shotgun sequencing, we assembled whole mitochondrial DNA genomes and estimated the timing of intraspecific divergence. Our projections of ENMs suggest current possum populations were unlikely to have been in contact during the Pleistocene. Although lowered sea level during glacial periods enabled connection with habitat in Tasmania, climate fluctuation during this time would not have facilitated gene flow over much of Australia. The most recent common ancestor of sampled intraspecific diversity dates to the early Pliocene when continental aridification caused significant changes to Australian ecology and Trichosurus vulpecula distribution was likely fragmented. Phylogenetic analysis revealed that the subspecies T. v. hypoleucus (koomal; southwest), T. v. arnhemensis (langkurr; north), and T. v. vulpecula (bilda; southeast) correspond to distinct mitochondrial lineages. Despite little phenotypic differentiation, Trichosurus vulpecula populations probably experienced little gene flow with one another since the Pliocene, supporting the recognition of several subspecies and explaining their adaptations to the regional plant assemblages on which they feed.


| BACKG ROU N D
Genetic signatures of range expansion, gene flow, and former barriers provide information that can help our understanding of past and future biological response to environmental change.
Phylogeographic studies of Australian animals reveal several different responses to past climatic oscillations. Some studies show that similar to the northern hemisphere, Pleistocene climate cycles were major drivers of species' phylogeographical patterns (Bryant & Fuller, 2014;Strasburg et al., 2007). But the majority of phylogeographic studies of Australian taxa suggest that older climate fluctuations explain current genetic structure (Kuch et al., 2005;Oliver et al., 2010Oliver et al., , 2017, especially among marsupials (Macqueen et al., 2010;Potter et al., 2012;Rowe et al., 2008). A rich fossil record of Australian megafauna allows calibration of changes in faunal composition (Hocknull et al., 2020;Miller et al., 2005;Roberts et al., 2001;Saltré et al., 2016;Turney et al., 2008) and it has been suggested that natural climate cycling during the Pleistocene had little impact on mammal assemblages (Prideaux et al., 2007). To determine whether current patterns of genetic structure within a single species date to events that occurred more than 2 million years ago, we examined the phylogeography of the Australian brushtail possum Trichosurus vulpecula.
Among Australia's endemic marsupial fauna, the brushtail possum Trichosurus vulpecula is one of the most widespread with a geographic range spanning the continent of nearly 7.7 million km 2 .
Populations of brushtail possum exist at both the east and west extremities of Australia, some 4000 km apart, and the species spans the temperate south to the tropical north ( Figure 1). Brushtail possums are arboreal herbivores, and across their geographic range, they interact with regional floras in ecologically distinct regions such as Jarrah and Kerri eucalyptus forest in Southwest Australia, tropical rainforests of Queensland, and temperate Tasmanian vegetation (Kerle, 1984). Despite this wide geographic and ecological range, the existence of numerous Aboriginal names (Abbott, 2012), and subspecies names for geographic races, no reliable diagnostic traits distinguishing brushtail possum populations are apparent from skull morphometrics, allozymes, or chromosomes (Kerle et al., 1991).

T A X O N O M Y C L A S S I F I C A T I O N
Population ecology F I G U R E 1 Sampling locations in Australia and New Zealand of Trichosurus vulpecula brushtail possums used for phylogenetic analyses. Spot colors coded by putative subspecies identification. The number of samples providing mitogenome sequences for molecular dating (in parentheses) and number of samples providing D-loop sequence [in square brackets]. In lighter green is the putative current distribution of brushtail possum in Australia (Kerle, 2002). Gray brushtail possum inset credit: Tony Jewell. The six subspecies of Trichosurus vulpecula are recognized based on geographic location, size, and fur color (How & Kerle, 1995), and show considerable physiological variation across Australia . Brushtail possum populations differ in their exposure to plant chemical defenses as they forage on plants known to contain compounds toxic to mammals including Erythrophleum, Acacia, Eucalyptus, and Gastrolabium, and as a result, have evolved different tolerances to toxins. For example, brushtail possums in Southwest Australia (koomal) have an LD50 160 times higher than their cousins (bilda) in East Australia (Twigg et al., 1996) when exposed to potent toxin fluoroacetate (known as 1080; Leong et al., 2017).
This emergence of physiological adaptations specific to regional flora demonstrates the coevolutionary interaction resulting in genomic divergence of spatial populations (Mead et al., 1985;Oliver & King, 1979), which implies a sustained interaction counter to the impression drawn from morphology.
Although Trichosurus vulpecula occurs widely in Australia, population density is uneven; there are large areas in the west and south of the country where the species is not recorded ( Figure 2a) and this could help explain lineage/subspecies distinctions. Lack of gene flow between brushtail possum populations in different habitats would allow fixing of local adaptations (Mallet, 2008). We predicted that past climate variation resulted in range shifts among brushtail possum populations across the Australian continent that resulted in a connection of habitat and populations that are now isolated.
Using modern records of brushtail possum occurrence, we infer the environmental envelope of Trichosurus vulpecula and extrapolate past distributions based on models of past climate. We then use independent time-calibrated phylogenetic analysis to test for the expected correlation between the timing of niche fragmentation and lineage splitting using a lineage-specific rate of molecular evolution based on a fossil-calibrated analysis.

| Niche modeling
To model the environmental envelope of Trichosurus vulpecula, we obtained occurrence data for the species in Australia from the Atlas of Living Australia (ala.org.au) which provided us with data from 80 different datasets (Appendix S1) including museum records from regions where they are not present anymore (dating as early as 1857; Kerle, 2002). Filtering duplicate records at a 0.1° latitude and longitude resolution resulted in a dataset of 6297 records.
F I G U R E 2 (a) Presence-absence points of Trichosurus vulpecula in Australia used for environmental niche modeling, with presence records from the GBIF database and pseudo-absence generated using the biomod2 SRE algorithm. Numbered squares represent the spatial blocks defined by blockCV and run separately to be combined in the ensemble modeling. (b) Weighted mean ensemble ecological niche modeling projections on current and past climate data (Mid-Holocene, Last Glacial Maximum, Last Interglacial period, Mid-Pliocene warm period, and Pliocene M2). Pixel of colors represents ensemble modeling scores; green, scores above the cut-off maximizing sensitivity (635.5); purple, scores >70% of this cut-off value (444.85); and gray, low scores signifying low probability of presence. Niche modeling approaches typically require both presence and absence data, but true absence information is difficult to verify.

Current
In this analysis, we used an alternative to true absence records by generating pseudo-absence data using the surface range envelop algorithm (SRE) from R package "Biomod2" (Thuiller et al., 2015). A total of 3000 pseudo-absence points were generated in locations and these were obtained from the Worldclim and Paleoclim databases (Brown et al., 2018;Fick & Hijmans, 2017). Multivariate environmental similarity surfaces analysis (MESS) was computed using the "dismo" R package  to check for the presence of non-analogous conditions in all scenarios of projection and to select only appropriate predictors that would not require models to extrapolate into novel climates. Multicollinearity of predictors was assessed using variance inflation factors (VIF, Zuur et al., 2010), retaining predictors with VIF smaller than 10 using "usdm" R package (Naimi et al., 2014). To investigate potential local differences and encompass Australia's broad climatic diversity and subspecies delimitations, the presence-absence data were separated into four spatial blocks using the R package "BlockCV" (Valavi et al., 2019;Figure 2a, T. v. fuliginosus and T. v. vulpecula were integrated in the same block, and T. v. johnsonii with T. v. eburacensis). The four spatial blocks were iteratively removed from the ecological niche models. We used a genetic algorithm from the R package "SDMtune" (Vignali et al., 2020) to optimize the models' hyperparameter configuration and obtain more robust models (Gradient Boosting Machine (GBM), Random Forest (RF), Artificial Neural Network (ANN), and Maximum Entropy (Maxent); Appendix S1: Table S1). The models were then computed using "Biomod2" in R (Thuiller et al., 2016)   , and added to this study in order to increase the sample size of the population from Tasmania and Victoria, which is their origin (Pracy, 1974). Total genomic DNA was extracted from 106 specimens using the GeneAid Tissue DNA Isolation Kit following the manufacturers' guidelines ( Figure 1), and quantity and quality checked using PerkinElmer LabChip® GX Touch HT.

| Acquisition of DNA sequence data
We examined the genealogy and population diversity of  Table S3).
D-loop sequences from the 17 Trichosurus vulpecula were extracted from the mitogenomes and included in the following analysis using this marker.
Discriminant analysis of principle components (DAPC: Jombart et al., 2010) was perform on the D-loop alignment using the "adegenet" package in R (Jombart, 2008) in order to describe the different population clusters. DAPC uses principal component analysis to select the components holding the most variance followed and then applies discriminant analysis that navigates among this variability to provide the clearest distinction between groups. We retained components that defined a total of 90% of the variance and we used the subspecies as the prior grouping factor for the analysis.

| Phylogenetic analysis
The newly assembled mtDNA genomes from 17 Trichosurus vulpecula were aligned with two mtDNA genomes from T. caninus. Phylogenetic relationships among these were inferred using maximum-likelihood and Bayesian inference methods using Raxml 8.2.12 (Stamatakis, 2014) and MrBayes v3.2.7 (Ronquist & Huelsenbeck, 2003). First, the data were divided into four partitions separating RNA (12S, 16S, and tRNA) and each codon position of protein-coding genes. PartitionFinder v2.1.1 (Lanfear et al., 2017) was then used to select the best-fitting substitution model for the data for each partition under corrected Akaike information criterion (AICc). This analysis determined the best substitution models to be GTR + G for all partitions and we used Raxml to compute the maximum-likelihood analyses with the default tree search approach using a simultaneous nearest-neighbor interchange method and a neighbor-joining tree as starting tree to estimate the ML tree topologies with 1000 bootstrap replicates (Felsenstein, 1985;Felsenstein & Kishino, 1993;Huelsenbeck & Hillis, 1993). Bayesian inference analyses were performed using unlinked parameters for each parti-

| Time calibration
We created an alignment of the mitochondrial genome from 14 species representing Dasyuromorphia and Diprotodontia, the two biggest Australian marsupials' orders (NCBI GenBank; Appendix S2: A molecular dating analysis was performed with this alignment to generate branch rate estimates that could be applied to the dating analysis of intraspecific Trichosurus lineages. The data were partitioned into tRNA and rRNA, D-loop, CDS first and second codon positions, and CDS third codon positions using Beauti  to input the priors and parameters for the Beast model. Substitution models for each partition were determined using the "BeastModel" algorithm. Model parameters included an optimized relaxed clock to allow the substitution rates to vary among the different species. A Yule model was implemented as it is appropriate for a phylogeny of mammal species as it assumes lineages split (constant speciation rate) without extinctions (Reid & Carstens, 2012;Steel & McKenzie, 2001). For time calibration of this phylogeny, we added priors estimated in a previous analyses of marsupial evolution using fossil calibrations (Meredith et al., 2009). That analysis of a DNA sequence from five nuclear genes (protein-coding portions of ApoB, BRCA1, IRBP, Rag1, and vWF) representing nine placental lineages and 22 marsupials included 32 hard time constraints based on both the fossil record and previous phylogenetic analyses (Meredith et al., 2008).
We applied mean times and 95% highest posterior density (HDP; to accommodate error around past estimates) to the nodes in our

For this, we implemented a coalescent constant population model in
Beast with an optimized relaxed clock to determine the mutation rate of each partition in the species', and a strict clock was then set to produce the final dated phylogeny. We used Treeannotator (Helfrich et al., 2019) to compile the trees from the Beast runs into a maximum clade credibility tree using a standard burn-in of the first 10% of trees and median node heights. Densitree and Figtree were used to visualize trees (Bouckaert, 2010).

| Ecological niche modeling
After filtering presence data, we had 5757 locality records of brushtail possums across Australia (Figure 2a). Mapping these records revealed five regions where living possums are documented in Australia, separated by areas of apparently uninhabited land (or sea). We retained five bioclimatic variables (following MESS and VIF analyses) that could be used over the region of study and over the four regions and six climatic conditions. Other bioclimatic variables were discarded to avoid extrapolating the models to never explored values and to improve robustness of models and projections (Appendix S1: Figure S1). Of the 64 runs, 55 had a ROC score greater than 0.80 and were retained for ensemble modeling (Appendix S1: Table S2). The ensemble model showed a good fit to the current range of Trichosurus vulpecula, with a cut-off value of 635.5 that maximized specificity (>98%) but still had a very good sensitivity of >93%. Plotting the projected ENM used this cut-off value (635.5) and 70% of this value to represent the available niche space (Figure 2b).

| Haplotypic diversity
Aligned D-loop (Appendix S2: Table S4) sequence data comprised 561 base pairs for 106 specimens and revealed 44 haplotypes differing from one another by up to 7% (uncorrected). The haplotype network revealed clusters broadly concordant with sampling locations ( Figure 3a). These clusters also appear in the DAPC analysis ( Figure 3b)

| Whole mitochondrial sequences
The mitochondrial genome of Trichosurus vulpecula was 14,701 bp, containing the expected 37 genes (13 protein-coding, 2 rRNA, and 22 tRNA). We aligned the 17 whole mtDNA sequences with two assembled from the related possum species Trichosurus caninus (Figure 1).
Our phylogenetic analyses resolved the same strongly supported topology (Appendix S2: Figure S4) with three well-supported lineages within our sampling of brushtail possums (maximum likelihood bootstrap score > 99 and Bayesian posterior probabilities >0.99).

Samples from Southwest Australia (Perth and Collie) and Barrow
Island that represent the subspecies T.

| Molecular dating
We calibrated a phylogenetic tree

| DISCUSS ION
In sexually reproducing animals, the fixation of adaptive traits in a population depends on the strength of natural selection and degree of gene flow (Mallet, 2008;Neher et al., 2010). While population size, range, and gene flow can be influenced by the emergence of adaptive traits (Reznick & Ghalambor, 2001;Szücs et al., 2017), climate and other abiotic factors have a strong influence on population structure (Davis & Shaw, 2001;Hewitt, 2004) and thus speciation (Endler, 1977). As a result, patterns of biodiversity and population structure are intrinsically linked to changes in Earth's climate (Byrne, 2008;Koot et al., 2022;Nogués-Bravo et al., 2018). Here, we tested the origin of intraspecific division in the Australian brushtail possum that experiences very different habitats across its spatial range. The signatures of Pliocene divergence have been observed in a number of widespread Australian animals (Macqueen et al., 2010;Oliver et al., 2017;Potter et al., 2012;Rowe et al., 2008). Here, we not only infer that Pliocene fragmentation has shaped the current brushtail possums' population structure but also suggest that Pleistocene glacial cycles have had an impact on their current distribution (Bryant & Fuller, 2014;Chapple et al., 2005).
Regional biota across the vast landscape of Australia have distinct characteristics such that Southwest Australia is a recognized biodiversity hotspot of global significance with endemic plants that account for at least 1.4% of the world's plant species (Myers et al., 2000;Rix et al., 2015), and forests of East Australia are also recognized as globally significant for their distinct biodiversity (Williams et al., 2011). Indeed, seven major ecoregions are recognized across the continent, each having distinct environmental attributes (Mittermeier, 2004), and a high level of endemism among plant and animal assemblages suggest protracted regional coevolution (Braithwaite, 1990;Firman et al., 2020;New, 2011;Porder, 2014). during earlier glacial periods (e.g., marine isotope stages 7 and 9), suggesting they had not diverged sufficiently to prevent gene flow.
Successful mixing of these two subspecies in New Zealand supports the idea that the two lineages have been isolated for brief periods followed by repeated connections maintaining reproductive compatibility (Pattabiraman et al., 2022). Perhaps the Tasmanian lineage could be better considered isolation by distance, rather than a "vicariant isolation" (Burridge, 2012).  during the Pliocene. Our molecular dating analysis using whole mitochondrial genome data suggests that the most recent common ancestor of the sampled populations of brushtail possum in Australia existed in the middle of the Pliocene (about 3.5 million years ago; 95% HDP interval 2.4, 4.6 Ma).
The phylogenetic results of this study would benefit from additional samples to include all six subspecies and samples from the Central and Eastern Australian populations. We have inferred phylogenies that are gene trees not subspecies trees, but there is no evidence that data from nuclear markers would change inferences about mtDNA lineage divergence. The concordant results from niche modeling and mtDNA phylogenetic analyses allow us to formulate an evolutionary scenario for brushtail possum in Australia ( Figure 5). The divergence between T. caninus and T. vulpecula corresponds to approximately the end of the Miocene, a period that is marked by a dry and cold climate in Australia (Zachos et al., 2001).
At the end of the Miocene, it is thought that rainforest persisted only in small patches in Southern and Southeastern Australia, with wet and open sclerophyll forests gradually replacing the rainforest (Black et al., 2012;Martin, 1998). This period of fragmented habitat might have been suitable for allopatric divergence of arboreal possums in their respective forest patches. Indeed, it appears that a F I G U R E 5 Hypothesis of intraspecific Trichosurus vulpecula differentiation in Australia inferred from calibrated phylogeny and niche modeling analyses. The main lineage splitting during the Pliocene corresponds to strong continental aridification, while Pleistocene glacial periods enabled exchange between Tasmania and South Australia due to lowered sea level. Aridification associated with the glacial period also restricted forest habitat that may have forced more intense plant-herbivore coevolution. Arid shrublands and grasslands emerge in northwestern central Australia.
Trichosurus vulpecula range retrenchment with regional populations increasingly isolated.
Climate cycling with pronounced fluctuations in precipitation. The overall trend throughout the Pleistocene was towards reduced precipitation.
Oscilating pattern of forest, grassland and desert habitat. Coastal land bridging during glacial periods, notably between Tasmania and south Australia.
Trichosurus vulpecula subject to local range restriction and expansion, including colonisation of Tasmania.

fu lg in o su s
Coevolution with regional flora.
number of modern marsupial groups (e.g., macropodine kangaroos, dasyurids, and peramelids) underwent extensive radiations during the Late Miocene, probably in response to this climate change (Beck, 2008;Meredith et al., 2009). Although late Miocene fossil deposits are rare in Australia, Trichosurus fossils have been found in the Riversleigh deposit (Archer et al., 2006;Roberts et al., 2008) confirming their presence in Eastern Australia at that time.
If the late Miocene can be considered a time of reduced niche space for Australian possums, then the early Pliocene provided more favorable conditions for brushtail possum population expansion. The short-lived Pliocene reversal of aridification yielded warmer, wetter conditions across most of Australia (Sniderman et al., 2016). Pollen records show that rainforests as well as wet sclerophyll forests expanded in peripheral regions (Black et al., 2012;Martin, 2006;Mcgowran & Bock, 2000) including semiarid Southern Australia (Sniderman et al., 2016). Development of these types of habitats at this time plausibly enabled expansion of T. vulpecula's range around the drier heartland of Australia, so that they reached the west from the south and/or the north. Fossils of a related species, Trichosurus hamiltonensis (Flannery et al., 1987) found near Hamilton, Victoria, indicate that suitable habitat for brushtail possums was available there at that time.
The Late Pliocene saw increasingly dry conditions but much of Collection CSIRO. We also thank OSPRI New Zealand and Predator Free New Zealand for the financing.

FU N D I N G I N FO R M ATI O N
Financial support for this research came from OSPRI New Zealand Ltd (Richard Curtis) and Predator Free 2050 Ltd (Dan Tompkins).

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare. All co-authors have seen and agree with the contents of the manuscript and there is no financial interest to report. We certify that the submission is original work and is not under review at any other publication.